#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Analysis and plot script for the output of "run/pde.py".
#

import os, sys, glob
import numpy as np
import matplotlib.pyplot as plt
import sqlite3

def wavelength(energy):
    '''
    @param energy - the energy in eV.
    @return wavelenght in nanometer.
    '''
    return 4.13566e-12 * 299.792458 / (energy * 1e-6)

files = glob.glob("./results/pde/*.sqlite")

wvl = []
pde = []
sigma_pde = []
for file in files:
    try:
        con = sqlite3.connect(file)
        cur = con.cursor()
        eMin, nParticles = cur.execute("SELECT eMin, nParticles FROM particleSourceMessenger;").fetchone()
        nHits = cur.execute("SELECT COUNT(*) FROM `g4sipmHits-0`;").fetchone()[0]
        wvl.append(wavelength(eMin))
        pde.append(nHits / float(nParticles) * 100.0)
        sigma_pde.append(np.sqrt(nHits) / float(nParticles) * 100.0)
    except:
        pass
       
# Sort results. 
idx = np.argsort(wvl)
wvl = np.asarray(wvl)[idx]
pde = np.asarray(pde)[idx]
sigma_pde = np.asarray(sigma_pde)[idx]

plt.plot(wvl, pde, '-', label="G4Sipm simulation")
plt.fill_between(wvl, pde - sigma_pde, pde + sigma_pde, color='#1f77b4', alpha=0.2)
plt.xlabel(r"wavelength $\lambda$ / nm")
plt.ylabel("photon detection efficicency / %")

# Rememeber axis limits.
xlim = plt.xlim()
ylim = plt.ylim()

# Reference values for the Hamamatsu S1036211100C from Hamamatsu
label = "Hamamatsu Datasheet (scaled by $c = 0.66$)"
wvl = [ 319.49, 320.882, 322.274, 323.666, 326.45, 327.842, 329.234, 330.626,
    332.019, 332.019, 333.411, 334.803, 336.195, 337.587, 338.979, 340.371, 343.155,
    344.548, 345.94, 348.724, 351.508, 355.684, 358.469, 364.037, 368.213, 372.39,
    375.174, 377.958, 382.135, 384.919, 387.703, 390.487, 393.271, 397.448, 400.232,
    404.408, 408.585, 412.761, 418.329, 422.506, 426.682, 432.251, 436.427, 441.995,
    447.564, 453.132, 458.701, 464.269, 468.445, 474.014, 479.582, 483.759, 487.935,
    493.503, 497.68, 501.856, 506.032, 510.209, 514.385, 518.561, 522.738, 526.914,
    529.698, 533.875, 538.051, 542.227, 546.404, 550.58, 553.364, 557.541, 561.717,
    565.893, 568.677, 572.854, 577.03, 581.206, 585.383, 588.167, 592.343, 596.52,
    600.696, 604.872, 609.049, 611.833, 616.009, 620.186, 624.362, 628.538, 632.715,
    636.891, 641.067, 645.244, 649.42, 653.596, 657.773, 661.949, 667.517, 671.694,
    675.87, 680.046, 684.223, 689.791, 693.968, 699.536, 703.712, 707.889, 713.457,
    717.633, 723.202, 727.378, 732.947, 737.123, 742.691, 746.868, 752.436, 756.613,
    762.181, 767.749, 771.926, 777.494, 781.671, 787.239, 792.807, 796.984, 802.552,
    808.121, 812.297, 817.865, 823.434, 827.61, 833.179, 838.747, 842.923, 848.492,
    854.06, 859.629, 865.197, 869.374, 874.942, 880.51, 886.079, 890.255, 895.824 ]

pde = [ 9.18188, 9.79401, 10.25308, 12.70162, 14.38495, 14.99709, 16.22136,
    16.83349, 18.05769, 18.05769, 19.28196, 19.89409, 20.9653, 21.57743, 22.18956,
    22.8017, 23.41383, 24.02597, 24.48503, 25.09717, 25.55624, 26.01537, 26.47444,
    26.93357, 27.23957, 27.69871, 28.15777, 28.61691, 29.07597, 29.68811, 30.14718,
    30.60631, 31.21844, 31.67751, 32.13658, 32.44265, 32.90178, 33.20785, 33.66691,
    33.97298, 34.27905, 34.58512, 34.89118, 35.04418, 35.04418, 35.04418, 34.89118,
    34.89118, 34.58512, 34.27905, 33.97298, 33.66691, 33.36085, 33.05478, 32.74871,
    32.28965, 31.98358, 31.52451, 31.21844, 30.75931, 30.30024, 29.84111, 29.53504,
    29.07597, 28.61691, 28.15777, 27.85171, 27.39264, 26.93357, 26.47444, 26.16837,
    25.7093, 25.25017, 24.7911, 24.33203, 23.8729, 23.41383, 23.10777, 22.64863,
    22.18956, 21.7305, 21.27136, 20.8123, 20.50623, 20.04709, 19.58803, 19.28196,
    18.82289, 18.36376, 18.05769, 17.59862, 17.13956, 16.83349, 16.37436, 16.06829,
    15.76222, 15.30315, 14.99709, 14.53802, 14.23195, 13.92589, 13.61982, 13.16069,
    13.00768, 12.54855, 12.24255, 11.93648, 11.63042, 11.32435, 11.01828, 10.71221,
    10.40615, 10.25308, 9.94701, 9.64101, 9.33495, 9.18188, 8.87581, 8.56975,
    8.26368, 7.95761, 7.80461, 7.49854, 7.19248, 7.03948, 6.73341, 6.58036, 6.27429,
    6.12126, 5.81519, 5.66216, 5.3561, 5.20307, 5.05004, 4.74397, 4.59095, 4.43791,
    4.13185, 3.97882, 3.82579, 3.67276, 3.36669, 3.21366 ]
plt.plot(wvl, pde, '.', label=label)

# Reference values for the Hamamatsu S10e6211100C from Tadday et. al.
label = "A. Tadday et. al. 2010"
wvl = [280, 350, 360, 370, 380, 390, 400, 410, 420, 430, 440, 450, 460, 470, 480, 490, 500,
            510, 520, 530, 540, 550, 560, 570, 580, 590, 600, 610, 620, 630, 640, 650, 660, 670, 680, 690, 700, 710,
            720, 730, 740, 750, 760, 770, 780, 790, 800, 810, 820, 830, 840, 850, 860, 870, 880, 890, 900, 910, 920,
            930, 940, 950, 960, 970, 980, 990, 1000]
pde = [0, 35, 33.25, 33, 33.5, 34.5, 34.95, 35.5, 35.6, 36.9, 36.9, 36.7, 35.95, 34.3, 35.1,
            34.5, 33.95, 33.05, 31.9, 30.75, 29.5, 28.4, 26.9, 25.6, 24.65, 23.55, 22.55, 21.55, 20.55, 19.5, 18.66,
            17.75, 16.9, 16.2, 15.33, 14.5, 13.8, 13.05, 12.3, 11.5, 10.9, 10.2, 9.55, 8.75, 8.15, 7.66, 7.1, 6.5, 6.05,
            5.66, 5.1, 4.6, 4.2, 3.85, 3.6, 3.2, 2.85, 2.5, 2.33, 2, 1.8, 1.6, 1.3, 1.05, 1, 0.75, 0.55,]
sigmaPde = [0, 2.5, 2.4, 2.4, 2.5, 2.5, 2.35, 2.4, 2.45, 2.5, 2.4, 2.3, 2.4, 2.25,
            2.3, 2.25, 2.1, 2.1, 2, 1.9, 1.8, 1.8, 1.7, 1.6, 1.55, 1.5, 1.45, 1.35, 1.35,
            1.25, 1.2, 1.05, 1.05, 1, 0.9, 0.9, 0.85, 0.8, 0.75, 0.75, 0.7, 0.7, 0.6, 0.55,
            0.5, 0.5, 0.4, 0.4, 0.4, 0.35, 0.35, 0.3, 0.3, 0.25, 0.25, 0.25, 0.25, 0.2, 0.2,
            0.15, 0.15, 0.1, 0.1, 0.05, 0.05, 0.02, 0.02 ]
# plt.errorbar(wvl, pde, fmt='o', yerr=sigmaPde, color="orange", label=label)

# Reference values for the Hamamatsu S10e6211100C from Tadday et. al.
label = "M. Lauscher, Aachen 2012"
wvl = [460, 402, 391, 375]
pde = [35.03200, 33.05468 , 30.86363, 26.86952]
sigmaPde = [3.3, 3.3, 3.3, 3.3, ]
plt.errorbar(wvl, pde, fmt='o', yerr=sigmaPde, label=label)

name, pitch, cellPitch, numberOfCells = cur.execute("SELECT name, pitch, cellPitch, numberOfCells FROM sipmModel;").fetchone()
plt.text(0.975, 0.975, "%s\n%d x %d mm, %d $\mu$m pitch\n%d cells"  % (name, pitch, pitch, cellPitch * 1000, numberOfCells), 
        ha="right", va="top", fontsize="medium", transform=plt.axes().transAxes)

plt.legend(loc="lower left")
plt.xlim(340, 640)
plt.ylim(18, 40)
plt.savefig("pde.pdf")
plt.show()
